Scalar self-force on eccentric geodesies in Schwarzschild spacetime: A time-domain 
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We calculate the self-force acting on a particle with scalar charge moving on a generic geodesic 
around a Schwarzschild black hole. This calculation requires an accurate computation of the retarded 
scalar field produced by the moving charge; this is done numerically with the help of a fourth-order 
convergent finite-difference scheme formulated in the time domain. The calculation also requires 
a regularization procedure, because the retarded field is singular on the particle's world line; this 
is handled mode-by-mode via the mode-sum regularization scheme first introduced by Barack and 
Ori. This paper presents the numerical method, various numerical tests, and a sample of results for 
mildly eccentric orbits as well as "zoom-whirl" orbits. 
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I. INTRODUCTION 

The inspiral and capture of solar-mass compact objects 
by supermassive black holes is one of the most promis- 
ing and interesting sources of gravitational radiation to 
be detected by the future space-based gravitational-wave 
antenna LISA [l[. For these extreme mass-ratio inspirals, 
one can treat the compact object as a point mass and de- 
scribe its influence on the spacetime perturbatively. Go- 
ing beyond the test mass limit, its motion is no longer 
along a geodesic of the unperturbed spacetime of the cen- 
tral black hole; it is a geodesic of the perturbed space- 
time created by the presence of the moving body. When 
viewed from the unperturbed spacetime, the small body 
is said to move under the influence of its gravitational 
self- force. The self- force induces radiative losses of energy 
and angular momentum, which will eventually drive the 
object into the black hole. To describe the motion of the 
body, including its inspiral toward the black hole, we seek 
to evaluate the self-force and calculate its effect on the 
motion. One way of doing this uses the mode-sum reg- 
ularization procedure introduced by Barack and Ori 
(For a comprehensive introduction of the problem, see 
the special issue of Classical and Quantum Gravity 0].) 

In this paper, in an effort to build expertise to calculate 
the gravitational self-force, we retreat to the technically 
simpler problem of a point particle of mass m endowed 
with a scalar charge q orbiting a Schwarzschild black hole 
of mass M. Following up on a previous paper [4|, we 
implement the numerical part of the regularization pro- 
cedure for generic orbits with a time-domain integration 
of the scalar-wave equation. 



A. The problem 

Our goal is to calculate the regularized self-force acting 
on a scalar point charge in orbit around a Schwarzschild 
black hole. In analogy with the gravitational case, where 
in a first-order (in m/M) perturbative calculation the 



particle moves on a geodesic of the background space- 
time, we take the orbit of the particle to be a geodesic and 
calculate the self-force as a vector field on this geodesic. 
We start by writing the Schwarzschild metric using the 
tortoise coordinate r* = r + 2Mln (^g — l) as 



ds 2 = / -dt 2 + dr* 



r 2 dfl 2 , 



(1.1) 



sm 2 ed(/) 2 ) is the 



where / = (l - 2M) , dtt 2 = (d 

metric on a two-sphere, and t, r, and <fi are the usual 
Schwarzschild coordinates. Our task is to solve the scalar 
wave equation 



g ap V a V p<$>(x) = -4^(x) 



fi(x) = q 5i(x,z(T))dT, 

J-y 



(1.2) 
(1.3) 



where V Q is the covariant derivative compatible with the 
metric g a fi, &(x) is the scalar field created by a scalar 
charge q which moves along a world line 7 : r 1— > z(t) 
parametrized by proper time r. The source term fx(x) 
appearing on the right-hand side is written in terms of a 
scalarized four-dimensional Dirac (5-function 8i{x,x') := 
S(x - Xq)5(xi - x' 1 )8(x 2 - x' 2 )S(x 3 - x' 3 )/*y-det(g a p). 
Because of the singularity in the source term, the re- 
tarded solution to Eq. (|1.2[) is singular on the world line, 
and the naive expression for the self-force, 



F a (r) = gV a *(z(r)), 



(1.4) 



must be regularized. Following DeWitt and Brehme [51, 
Mino, Sasaki, Tanaka 0, Quinn and Wald Q, Quinn [g] 
carried out this regularization for the electromagnetic, 
scalar and gravitational radiation reaction. In later work, 
Detweiler and Whiting [9( introduced a very useful de- 
composition of the retarded solution of Eq. (|1.2p in terms 
of a singular part $> s and a regular remainder & R : 



(1.5) 



$ fl is regular and differentiable at the position of the par- 
ticle, satisfies the homogeneous wave equation associated 
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with Eq. (|1.2p , and is solely responsible for the self- force 
acting on the particle. $ s , on the other hand, satisfies 
Eq. (| 1 . 2(1 . is just as singular at the particle's position as 
the retarded solution, and produces no force on the par- 
ticle. Rearranging Eq. (|1.5p and differentiating once, we 
can write the regularized self-force as 



(1.6) 



In a previous paper jj], we described our implemen- 
tation of the regularization procedure to find a mode- 
sum representation of V a $ along a generic geodesic of 
the Schwarzschild spacetime. Schematically, we intro- 
duce a tetrad e? •> and decompose the tetrad components 
:= e?iV a $ of the field gradient in terms of ordinary 
scalar spherical harmonics Y* 



( III • 



(1.7) 



Each mode $^j(t, r) is finite at the position of the par- 
ticle, but their sum diverges on the world line. In we 
derive analytic expressions for the mode-sum decompo- 
sition of 
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C, 



l+h 



D 



(1.8) 



(1.9) 



where the coefficients Ar^ , B(jA , CV^) , and Di^ are in- 
dependent of £; they are listed in Appendix [B] for conve- 
nience. 

As each mode of $ is finite, it is straightforward to 
compute the modes of the retarded solution using nu- 
merical methods, and we will describe how this was done 
in Sec. IIV1 We use the numerical solutions in Eq. (|1.6p 
to calculate the regularized self-force, regularizing mode- 
by-mode: 



00 



E< 



)> 



(1.10) 



where < I > ( AI ).£ := ^2 m ^uK^tm ( n0 summation over £ im- 
plied). 

For numerical purposes it is convenient to define ipi m 
by 



00 t 1 

$ W = E E -^m^™, 



(1.11) 



where Yi m are the usual scalar spherical harmonics. Af- 
ter substituting in Eq. (|1.2p , this yields a reduced wave 
equation for the multipole moments ipe m - 

-dflpim + c%»ip im - Vnptm = 

" 4:Trq^=Y im (7r/2, cj> )6(r* - r*), (1.12) 



where 



V t = f 



/2M £(£+1) 



(1.13) 



An overbar denotes complex conjugation, _E = —ut is the 
particle's conserved energy per unit mass, and u a — ^— 
is its four velocity. Quantities bearing a subscript "0" are 
evaluated at the particle's position; they are functions of 
r that are obtained by solving the geodesic equation 



Van" = 



(1.14) 



in the background spacetime. Without loss of general- 
ity, we have confined the motion of the particle to the 
equatorial plane 9 = | . 

Once we have numerically solved Eq. (|1.12[) . we ex- 
tract numerical estimates for if)e m , dtipim and d r *?pe m , 
which can then be used to find $£ m , dt$i m and d r &£ m . 
These — together with the translation table displayed in 
Eqs. (1.23)-(1.26) of 0], reproduced in Appendix 0— 



allow us to find the tetrad components $ 



with re- 



spect to the tetrad defined by Eqs. (1.18)-(l.21) of [4j. 
Eventually we regularize the multipole coefficients 

i 

®W = E ^Me m (t ,ro)Y em (Tr/2 7 M (1.15) 



using Eq. (|1.10[) ; this involves the regularization param- 
eters listed in Eqs. (1.30)-(1.45) of (4j, which are repro- 
duced in Appendix [Bj 



B. Organization of this paper 

In Sec. |TT] we introduce the main ideas behind the 
discretization scheme used in the numerical simulation. 
Sec. Mil describes the choices we make in order to handle 
the problems of specifying initial data and proper bound- 
ary conditions. The next section — Sec. lIVI — provides de- 
tails on the concrete implementation of the ideas put 
forth in Sees. |TT]and|TTTJ In Sec. |V]we describe the tests 
we performed in order to validate our implementation of 
the numerical method. Sec. IVII finally presents sample 
results for a small number of representative simulations. 



C. Future work 

This work, which deals with a scalar charge moving in 
the Schwarzschild spacetime, is not intended to produce 
physically or astrophysically interesting results. Instead, 
its goal is to help us evaluate the merits of several strate- 
gies that could be used to tackle the more interesting (and 
difficult) problems of electromagnetism and gravity. 

One future project we are currently exploring is to ap- 
ply the formalism developed so far to the electromagnetic 
self-force acting on an electric charge. Beyond the tech- 
nical complication of having to deal with a vector field 
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instead of a single scalar quantity, we are also faced with 
the reality of having to impose a gauge (in our case: the 
Lorenz gauge) and to eliminate (or at least control) gauge 
violations in the numerical simulation. The first step, 
namely, the calculation of the regularization parameters 
AifA, B(jA, C(^), and for the self-force, is currently 
underway. Also underway is the calculation of the regu- 
larization parameters for he gravitational self-force. 

Another project is the implementation of a scheme to 
use the calculated self-force to update the orbital pa- 
rameters of a particle on its inspiral toward the black 
hole. The standard proposed approach to this problem 
in the past has been to calculate the self-force on a set 
of geodesies which are momentarily tangent to the par- 
ticle's trajectory. The self-force calculated in this way 
is then used to update the orbital elements. This "after 
the fact" calculation of the motion requires one to build 
(in advance) a large database of self-force values for the 
anticipated set of orbital parameters that the particle's 
trajectory will assume during its inspiral. Alternatively, 
and conceptually more simply, the self-force could be cal- 
culated self-consistently along the real, accelerated tra- 
jectory. Such an approach requires changes in the expres- 
sions of the regularization parameters, which so far have 
been derived only for geodesic orbits. We are currently 
investigating the merits of such an approach. 



II. NUMERICAL METHOD 

In this section we describe the algorithm used to inte- 
grate the reduced wave equation [Eq. (|1.12p ] numerically. 
For the most part we use the fourth-order algorithm in- 
troduced by Lousto [1(3] , with some modifications to suit 
our needs. We choose to implement a fourth-order con- 
vergent code because second-order convergence for the 
potential while much easier to achieve, would guaran- 
tee only first-order convergence for V Q $, the quantity in 
which we are ultimately interested. With a fourth-order 
convergent code we can expect to achieve third-order con- 
vergence for V Q <!>, which is required for an accurate es- 
timation of the self-force. Numerical experiments, how- 
ever, show that in practice we do achieve fourth-order 
convergence for the derivatives of a fortunate outcome 
that we exploit but cannot explain. 

From now on, we will suppress the subscripts i and m 
on Vi and ipt m for convenience of notation. The wave 
equation consists of three parts: the wave-operator term 
(d%» — df)ip and the potential term Vip on the left-hand 
side, and the source term on the right-hand side of the 
equation. Of these, the wave operator turns out to be 
easiest to handle, and the source term does not create a 
substantial difficulty The term involving the potential 
V turns out to be the most difficult one to handle. 

Following Lousto we introduce a staggered grid with 
step sizes At = \Ar* = h, which follows the characteris- 
tic lines of the wave operator in Schwarzschild spacetime; 
see Fig. [T] for a sketch of a typical grid cell. The basic 



idea behind the method is to integrate the wave equation 
over a unit cell of the grid, which nicely deals with the 
Dirac-<5 source term on the right-hand side. To this end, 
we introduce the Eddington-Finkclstcin null coordinates 
v = t + r* and u = t — r* and use them as integration 
variables. 



A. Differential operator 

Rewriting the wave operator in terms of u and v, we 
find — &l + d^, = —4d u d v , which allows us to evaluate 
the integral involving the wave operator exactly. We find 



cell 



-4d u d v ip dudv = - 4[ip(t + h,r*) + tp(t - h,r*) 

-ip(t,r* -h)-ip(t,r* + h)}. 

(2.1) 



B. Source term 

If we integrate over a cell traversed by the particle, then 
the source term on the right-hand side of the equation 
will have a non-zero contribution. Writing the source 
term as G(t,r*)S(r* - r^{t)) with 



G(t,r*) = -4Trq-t-Y irn {n/2,<j> ), 



(2.2) 



we find 



GS(r*-r *(*))dud W = -^r« 
ai & Jt t r o(t) 

x y €m (7r/2,^o(t))dt, 

(2.3) 

where t\ and t<i are the times at which the particle enters 
and leaves the cell, respectively. While we do not have 
an analytic expression for the trajectory of the particle 
(except when the particle follows a circular orbit), we can 
numerically integrate the first-order ordinary differential 
equations that govern the particle's motion to a precision 
that is much higher than that of the partial differential 
equation governing ip. In this sense we treat the integral 
over the source term as exact. To evaluate the integral 
we adopt a four-point Gauss-Legendre scheme, which has 
an error of order h s . 



C. Potential term 

The most problematic term — from the point of view of 
implementing an approximation of sufficiently high or- 
der in h — turns out to be the term Vip in Eq. p. 1 2|l . 
Since this term does not contain a 5-function, we have to 
approximate the double integral 
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Vip du dv 



(2.4) 
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FIG. 1: Points used to calculate the integral over the potential 
term for vacuum cells. Grid points are indicated by blue cir- 
cles while red cross-hairs indicate points in between two grid 
points. We calculate field values at points that do not lie on 
the _grid by employing the second-order algorithm described 
in Qil- 



up to terms of order ft 6 for a generic cell in order to 
achieve an overall 0(ft 4 ) convergence of the scheme. 

Here we have to treat cells traversed by the particle 
("sourced" cells) differently from the generic ("vacuum") 
cells. While much of the algorithm can be transferred 
from the vacuum cells to the sourced cells, some modifica- 
tions are required. We will describe each case separately 
in the following subsections. 



1. Vacuum case 

To implement Lousto's algorithm to evolve the field 
across the vacuum cells, we use a double Simpson rule to 
compute the integral Eq. (|2.4[) . We introduce the nota- 
tion 



g(t,r*) = V(r*)iP(t,r*) 



(2.5) 



and label our points in the same manner (see Fig. [T]) as 
in 0: 



gdudv=[~) [gi 

cell V 



52+33+54 
324 + 334 + 313) + 163o] H 



+■ 4(312+ 

0(ft 6 ). (2.6) 



Here, for example, 31 is the value of 3 at the grid point 
labeled 1, and 312 is the value of 3 at the off-grid point 
labeled 12, etc. Deviating from Lousto's algorithm, we 
choose to calculate 30 using an expression different from 
that derived in [io| . Unlike Lousto's approach, our ex- 
pression exclusively involves points that are within the 
past light cone of the current cell. We find 



3o 



1G 



[8V 4 V>4 + 8V1 A + 8V 2 ^ 2 - 4U 6 A - 4U 5 ^5 
Vw iho + V 7 i> 7 - Vq ^9 - V s ^s] + 0(ft 4 ). (2.7) 



In order to evaluate the term in parentheses in 
Eg. (|2.6I). we again use a variant of the equations given 
in Lousto's equations (33) and (34), 



3i3 + 312 =V(r* Q - ft/2) (fa + tpo) 



1- 1(| ) IV, 1 , -ft/2) 



324 + 334 =V(r* + ft/2) (V> + ipi) 



0(ft 4 ), 



(2.8) 



0(ft 4 ) 



(2.9) 



contain isolated occurrences of -i/ , o ) the value of the field 
at the central point. Since Eq. (|2.7[) only allows us to 
find 30 = VbV'Oi finding ipQ would involve a division by 
Vo, which will be numerically unstable very close to the 
event horizon where Vo ~ 0. Instead we choose to express 
the potential term appearing in the square brackets as a 
Taylor series around 7q. This allows us to eliminate the 
isolated occurrences of ipo, and we find 



313+312 + 324 + 334 = 2V{r* )ip 
+ V(r*-h/2)^ 



1-5(5! ™> 



+ V(r* + ft/2) ^4 



+ \[V(r*-h/2)-2V(r* ) 

+ V(r* +h/2)] (V>!+^ 4 ) + 0(ft 4 ). 



(2.10) 



Because of the (3 



factor in Eq. (|2.6|) , this allows us 
to reach the required 0(ft 6 ) convergence for a generic 
vacuum cell. This — given that there is a number of order 
N = 1/ft 2 of such cells — yields the desired overall 0(ft 4 ) 
convergence of the full algorithm, at the end of the N 
steps required to finish the simulation. 



2. Sourced cells 

For vacuum cells, the algorithm described above is 
the complete algorithm used to evolve the field forward 
in time. For cells traversed by the particle, however, 
we have to reconsider the assumptions used in deriving 
Eqs. (j2~?) and (|2~TU|) . When deriving Eq. (|2~TU)l we have 
employed the second-order evolution algorithm (see [Ic| ) , 
in which the single step equation 



ip3 



1p2 



Vo (ipi + 4>a) 



(2.11) 
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FIG. 2: Cells affected by the passage of the particle, showing 
the reduced order of the single step equation 



is accurate only to 0(h 3 ) for cells traversed by the parti- 
cle. For these cells, therefore, the error term in Eq. (|2.10p 
is 0(h 3 ) instead of 0(h 4 ). As there is a number of or- 
der N' = 1/h of cells that are traversed by the parti- 
cle in a simulation run, the overall error — after including 
the (|) 2 factor in Eq. (|2.6p — is of order h 4 . We can 
therefore afford this reduction of the convergence order 
in Eq. (|2~TU|) 

Equation (|2.7|) . however, is accurate only to O(h) for 

cells traversed by the particle. Again taking the (-|) 
factor into account, this renders the overall algorithm 
0(h 2 ). Figure [2] shows the cells affected by the particle's 
traversal and the reduced order of the single step equa- 
tion for each cell. Cells whose convergence order is 0(h 5 ) 
or higher do not need modifications, since there is only a 
number N' = 1/h of such cells in the simulation. We are 
therefore concerned about cells neighboring the particle's 
trajectory and those traversed by the particle. 

a. Cells neighboring the particle These cells are not 
traversed by the particle, but the particle might have 
traversed cells in their past light-cone, which are used in 
the calculation of go in Eq. (12. 7| . For these cells, we use 
a one-dimensional Taylor expansion of g(t, r*) within the 
current time-slice t — t , 
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[5V(r* -h)i>(t ,r* -h) 

+ 15V(r*-3h)ij{t ,r*~3h) 

-5V(r*-5h)yj(t ,r*-5h) 

+ V(r* -7h)tP(t ,r*~7h)]+O(h 4 ) 

for the cell on the left-hand side, and 



(2.12) 



9o 



16 



[5V(r*+h)^(t ,r* + h) 
15V(r^ + 3h)ip(t ,r^ + 3h) 



-5V{r*+5h)iP{t ,r*+5h) 

+ V{r* + 7h) <tp(t , r* + 7h)] + 0(h 4 ) (2.13) 

for the cell on the right-hand side, where (to,^) ' s the 
center of the cell traversed by the particle. Both of these 
are more accurate than is strictly necessary; we would 
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FIG. 3: Typical cell traversal of the particle. We split the 
domain into sub-parts indicated by the dotted line based on 
the time the particle enters (at ti) and leaves (at fe) the cell. 
The integral over each sub-part is evaluated using an iterated 
two-by-two point Gauss-Legendre rule. 

need error terms of order h 3 to achieve the desired over- 
all 0(h 4 ) convergence of the algorithm. Keeping the ex- 
tra terms, however, improves the numerical convergence 
slightly. 

b. Cell traversed by the particle We choose not to 
implement a fully explicit algorithm to handle cells tra- 
versed by the particle, because this would increase the 
complexity of the algorithm by a significant factor. In- 
stead we use an iterative approach to evolve the field 
using the integrated wave equation 



-4(^3+^2 -ipi- ipi) 



V ip du dv = 



cell 



--ft/ ^Y em (ir/2,Mt))dt. (2.14) 



E J tl r Q (t) 

In this equation the integral involving the source term 
can be evaluated to any desired accuracy at the begin- 
ning of the iteration, because the motion of the particle 
is determined by a simple system of ordinary differential 
equations, which are easily integrated with reliable nu- 
merical methods. It remains to evaluate the integral over 
the potential term, which we do iteratively. Schemati- 
cally the method works as follows: 

• Make an initial guess for ip^ using the second-order 
algorithm. This guess is correct up to terms of 
0(h 3 ). 

• Match a second-order piecewise interpolation poly- 
nomial to the six points that make up the past light- 
cone of the future grid point, including the future 
point itself. 

• Use this approximation for tp to numerically calcu- 
late 



V tpdu dv, 



using two-by-two point Gauss-Legendre rules for 
the six sub-parts indicated in Fig. [3] 

• Update the future value of the field and repeat the 
process until the iteration has converged to a re- 
quired degree of accuracy. 
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FIG. 4: Numerical domain evolved during the simulation. We 
impose an inner boundary condition close to the black whole 
where we can implement it easily to the accuracy of the un- 
derlying floating point format. Far away from the black hole, 
we evolve the full domain of dependence of the initial data 
domain without imposing boundary conditions. 



III. INITIAL VALUES AND BOUNDARY 
CONDITIONS 

As is typical for numerical simulations, we have to pay 
careful attention to specifying initial data and appropri- 
ate boundary conditions. These aspects of the numerical 
method are highly non-trivial problems in full numerical 
relativity, but they can be solved or circumvented with 
moderate effort in the present work. 



outgoing boundary conditions at spatial infinity r* — > oo, 
ie. 



lim d u ip ~0, 



lim d v ip =0. 



(3.2) 



Because of the finite resources available to a computer 
we can only simulate a finite region of the spacetime, 
and are faced with the reality of implementing boundary 
conditions at finite values of r*. Two solutions to this 
problem present themselves: 

1. choose the numerical domain to be the domain of 
dependence of the initial data surface. Since the 
effect of the boundary condition can only propagate 
forward in time with at most the speed of light, 
this effectively hides any influence of the boundary. 
This is what we choose to do in order to deal with 
the outer boundary condition. 

2. implement boundary conditions sufficiently "far 
out" so that numerically there is no difference be- 
tween imposing the boundary condition there or at 
infinity. Since the boundary conditions depend on 
the vanishing of the potential V(r) appearing in the 
wave equation, this will happen once 1 — 2M / r»0. 
Near the horizon r w 2M(1 + exp(r*/2M)), so 
this will happen — to numerical accuracy — for mod- 
estly large (negative) values of r* — 73 M. We 
choose to implement the ingoing waves condition 
d u tpim = there. 



A. Initial data 

In this work we use a characteristic grid consisting of 
points lying on characteristic lines of the wave operator 
to evolve i\) forward in time. As such, we need to specify 
characteristic initial data on the lines u = uq and v = vq 
shown in Fig. 2] We choose not to worry about specifying 
"correct" initial data, but instead arbitrarily choose tp to 
vanish on u — uq and v — vq: 



■0(u = Uq, v) = l/j(u, V = Vq) = 0. 



(3.1) 



This is equivalent to adding spurious initial waves in the 
form of a homogeneous solution of Eq. (|1.12p to the cor- 
rect solution. This produces an initial wave burst that 
moves away from the particle with the speed of light, 
and quickly leaves the numerical domain. Any remain- 
ing tails of the spurious initial data decay as i _ ( 2£ + 2 ) as 
shown in and become negligible after a short time. 
We conclude that the influence of the initial-wave con- 
tent on the self-force becomes negligible after a time of 
the order of the light-crossing time of the particle's orbit. 



B. Boundary conditions 

On the analytical side we would like to impose ingoing 
boundary conditions at the event horizon r* — > — oo and 



IV. IMPLEMENTATION 

Making more precise the ideas developed in the pre- 
ceding sections, we implement the following numerical 
scheme. 



A. Particle motion 

Following Darwin [l2j we introduce the dimensionless 
semi-latus rectum p and the eccentricity e such that for 
a bound orbit around a Schwarzschild black hole of mass 
M, 



P M 
1 + e 



'2 



pM 
1 - e 



(4.1) 



are the radial positions of the periastron and apastron, 
respectively. Energy per unit mass and angular momen- 
tum per unit mass are then given by 



E 2 



{p-2- 2e)(p - 2 + 2e) 
p {p — 3 — e 2 ) 



p 2 M 2 



p — 3 — e 2 



(4.2) 



Together with these definitions it is useful to introduce 
an orbital parameter \ such that along the trajectory of 
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the particle, 



r{x) 



pM 



1 + e cos x 



(4.3) 



where \ is single- valued along the orbit. We can then 
write down first-order differential equations for x(t) and 
the azimuthal angle (f>(t) of the particle, 

dx (p — 2 — 2ecosx)(l + ecosx)(l + ecos%) 
~dt = 



(Mp 2 ) 



p — 6 — 2e cos x 



y (p-2-2e)(p-2 + 2e)' 

d0 (p — 2 — 2ecos%)(l + ecosx) 2 
dt ~ p 3/2M^/(p - 2 - 2e) (p - 2 + 2e) ' 



(4.4) 
(4.5) 



We use the embedded Runge-Kutta-Fehlberg (4, 5) algo- 
rithm provided by the GNU Sc i ent ific Library routine 
gsl_odeiv_stepjrkf 45 and an adaptive step-size control 
to evolve the position of the particle forward in time. 
Intermediate values of the particle's position are found 
using a Hermite interpolation of the nearest available cal- 
culated positions. 



B. Initial data 

We do not specify initial data. The field is set to zero 
on the initial characteristic slices, u = uq and v = Vq. 



C. Boundary conditions 

We adjust the outer boundary of the numerical do- 
main at each time-step so that we cover the domain of 
dependence of the initial characteristic surfaces and the 
particle's world line. The resulting numerical domain was 
already shown in Fig. [U 

Near the event horizon, at r* sa —73 M, we implement 
an ingoing-wave boundary condition by imposing 



iP(t + h,r*) = ip(t,r* -h). 



(4.6) 



This allows us to drastically reduce the number of cells 
in the numerical domain, and consequently the running 
time of the simulation. 



+ 



1 -l(l)Vo + Vi) + lQ)Vo(V„ + % ) 



^4 



1 (h 



4 V 3 



V 



(912 + 524 + 534 + 513 + 45o), 

(4.7) 



where go is given by Eq. (|2.7|) and the sum 1712 + 524 
534 + 5i3 is given by Eq. (|2. 10|) . 



E. Cells next to the particle 

Vacuum cells close to the current position of the parti- 
cle require a different approach to calculate 50, since the 
cells in their past light cone could have been traversed 
by the particle. We use Eqs. (j2~T2|) and (|2"T3"|l to find 
go in this case. Other than this modification, the same 
algorithm as for generic vacuum cells is used. 



F. Cells traversed by the particle 

We evolve cells traversed by the particle using the it- 
erative algorithm described in Sec. Ill C 21 Here 



+ 04 



■01 + i>2 
1 

4 



V tp du dv 



cell 



1p3 



E J H r Q (t) 

where the initial guess for the iterative evolution of 
J J Vip du dv is obtained using the second order algo- 
rithm of Lousto and Price [131 , 



r^y< m (7r/2,7ro(t))dt, (4.8) 



f/'3 



X [lp 2 + 'tjji] 
2nq [*> fa® 
E J tl r (t) 



Y im (n/2,n (t))dt. (4.9) 



Successive iterations use a four-point Gauss-Legendre 
rule to evaluate the integral of Vip; this requires a second- 
order polynomial interpolation of the current field values 
as described in Appendix [Cl 



D. Evolution in vacuum 



G. Extraction of the field data at the particle 



Cells not traversed by the particle are evolved using 
Eqs. (|2Tj) , (f276|) - IpTTU)) . Explicitly written out, we use 



1p3 = -1p2 



In order to extract the value of the field and its first 
derivatives at the position of the particle, we again use 
a polynomial interpolation at the points surrounding the 
particle's position. Using a fourth-order polynomial, as 
described in Appendix [Cl we can estimate ip, dtipt, and 
d r * ip at the position of the particle up to errors of order 
h A . As was briefly mentioned in Sec.HH we would expect 
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an error term of order h 3 for dtipt and d r *ip. The 0(h 4 ) 
accuracy we actually achieve by using a fourth-order (in- 
stead of a third-order) piecewise polynomial shows up 
clearly in a regression plot such as Fig. [7] 

H. Regularization of the mode sum 

We use the calculated multipole moments iptm to con- 
struct the multipole moments $£ m , and first derivatives 
dt^im and d r $im,, of the scalar field. These, in turn, are 
used to calculate the tetrad components $(o)i m , 3>(+)f m > 
$(_)i m , and $(3)t m of the field gradient according to 
Eqs. (1.23)-(1.26) of [1], which are reproduced in Ap- 
pendix [X] These multipoles then give rise to the multi- 
pole coefficients of the retarded field, 

t 

* M *(t,r,M)= *00*m(*» r )*WM), (4-10) 

which are subjected to the regularization procedure de- 
scribed by Eq. (1.29) of [|, 

$f M) (t,r o ,7r/2,0 o ) = hm V{$(^(*»^o+A,7r/2,^o) 
-q[(£ + l/2)A M + B w 

{£ + 1/2) + (£-l/2)(e + 3/2) 
+ ...]}, (4.11) 

using the regularization parameters ArjA, B^, an( i 
Z?( M ) tabulated in Appendix IB1 

Finally we reconstruct the vector components of the 
field gradient using Eqs. (1.47)-(1.48) of [|, 

$? = VTo$to), (4-12) 

= -= ($f +) e-^° + $f_ )e l< M , (4.13) 
v Jo ^ ' 

^ = -r $f 3) , (4.14) 
d>« = ($f +) e-^° - $f_ )e ^) , (4.15) 

and calculate the self-force 

F a = gtf* (4.16) 

We recall the discussion in Sec. II Al concerning the def- 
inition of $> R , its connection to the self- force acting on 
the particle, and its regularity at the particle's position. 

V. NUMERICAL TESTS 

In this section we present the tests we have performed 
to validate our numerical evolution code. First, in order 



to check the fourth-order convergence rate of the code, 
we perform regression runs with increasing resolution for 
both a vacuum test case, where we seeded the evolution 
with a Gaussian wave packet, and a case where a particle 
is present. As a second test, we compute the regularized 
self-force for several different combinations of orbital el- 
ements p and e and check that the multipole coefficients 
decay with I as expected. This provides a very sensi- 
tive check on the overall implementation of the numerical 
scheme, as well as the analytical calculations that lead to 
the regularization parameters. Finally, we calculate the 
self-force for a particle on a circular orbit and show that 
it agrees with the results presented in 0, [U| • 

A. Convergence tests: Vacuum 

As a first test of the validity of our numerical code we 
estimate the convergence order by removing the particle 
and performing regression runs for several resolutions. 
We use a Gaussian wave packet as initial data, 

ip{u = u ,v) =exp(-[v-u p ] 2 /[2cr 2 ]), (5.1) 
i>(u,v = v ) = 0, (5.2) 

where v p = 75 M and a — 10 M, vq — —uq = 6M + 
2 M In 2, and we extract the field values at r* = 20 M. 
Several such runs were performed, with varying resolu- 
tion of 2, 4, 8, 16, and 32 grid points per M. Figure [5] 
shows ip(2h) — ijj(h) rescaled by appropriate powers of 2, 
so that in the case of fourth-order convergence the curves 
would lie on top of each other. As can be seen from the 
plots, they do, and the vacuum portion of the code is 
indeed fourth-order convergent. 

B. Convergence tests: Particle 

While the convergence test described in section IV Al 
clearly shows that the desired convergence is achieved 
for vacuum evolution, it does not test the parts of the 
code that are used in the integration of the inhomoge- 
neous wave equation. To test these we perform a second 
set of regression runs, this time using a non-zero charge 
q. We extract the field at the position of the particle, 
thus also testing the implementation of the extraction 
algorithm described in section IIV Gl For this test we 
choose the I = 6, m = 4 mode of the field generated 
by a particle on a mildly eccentric geodesic orbit with 
p = 7, e = 0.3. As shown in Fig. [5] the convergence 
is still of fourth order, but the two curves no longer lie 
precisely on top of each other at all times. The region 
before t » 100 M is dominated by the initial wave burst 
and therefore does not scale as expected, yielding two 
very different curves. In the region 300 M < t < 400 M 
the two curves lie on top of each other, as expected for a 
fourth-order convergent algorithm. In the region between 
t w 200 M and t w 300 AT, however, the dashed curves 
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FIG. 5: Convergence test of the numerical algorithm in the 
vacuum case. We show differences between simulations using 
different step sizes h = 0.5 M (i/> 2 ), h = 0.25 M (<0 4 ), h = 
0.125 M (ip 8 ), h = 0.0625 M (Vie), and ft = 0.03125 M (V32). 
Displayed are the rescaled differences 84-2 — ip4 — ipi, 8^-4 = 
2 4 (V>8 - 1P4), <5i6-8 = 4 4 (i/;8 - V'4), and 832-16 = - ^4) 

for the real part of the £ = 2, m = 2 mode at r* ~ 20 M. 
The maximum value of the field itself is of the order of 0.1, 
so that the errors in the field values are roughly five orders 
of magnitude smaller than the field values themselves. We 
can see that the convergence is in fact of fourth-order, as the 
curves lie nearly on top of each other, with only the lowest 
resolution curve 84,-2 deviating slightly. 



have slightly smaller amplitudes than the solid one, indi- 
cating an order of convergence different from (but close 
to) four. 

To explain this behavior we have to examine the terms 
that contribute significantly to the error in the simula- 
tion. The numerical error is almost completely domi- 
nated by that of the approximation of the potential term 
//ceil ^V'dwdv in the integrated wave equation. For vac- 
uum cells the error in this approximation scales as h e , 
where h is the step size. For cells traversed by the parti- 
cle, on the other hand, the approximation error depends 
also on the difference t^— 1\ of the times at which the par- 
ticle enters and leaves the cell. This difference is bounded 
by h but does not necessarily scale as h. For example, if 
a particle enters a cell at its very left, then scaling h by \ 
would not change ti — t\ at all, thus leading to a scaling 
behavior that differs from expectation. 

To investigate this further we conducted test runs 
of the simulation for a particle on a circular orbit at 
r — 6 M. In order to observe the expected scaling behav- 
ior, we have to make sure that the particle passes through 
the tips of the cell it traverses. When this is the case, then 
t2 — t\ = h and a plot similar to the one shown in Fig. [6] 
shows the proper scaling behavior. As a further test we 
artificially reduced the convergence order of the vacuum 
algorithm to two by implementing the second-order algo- 
rithm described in [10|. By keeping the algorithm that 
deals with sourced cells unchanged, we reduced the rela- 



FIG. 6: Convergence test of the numerical algorithm in the 
sourced case. We show differences between simulations using 
different step sizes of 4 (t/>4), 8 (V's), 16 (i/>i6), and 32 (^32) 
cells per M. Displayed are the rescaled differences 8g,-4 = 
t/>8 — *p4, etc. (see caption of Fig. [S] for definitions) of the field 
values at the position of the particle for a simulation with 
£ — 6, m — 4 and p — 7, e = 0.3. We see that the convergence 
is approximately fourth-order. 
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FIG. 7: Convergence test of the numerical algorithm in the 
sourced case. We show differences between d r & for simula- 
tions using different step sizes of 4 ($r-, 4 ), 8 (^V.s), 16 {& r ,i6), 
and 32 ($ r ,32) cells per M. Displayed are the rescaled differ- 
ences 8s-4 = <&r,8 — 3v,4 etc. of the values at the position of 
the particle for a simulation with £ — 6, m = 4 and p = 7, 
e = 0.3. Although there is much noise caused by the piece- 
wise polynomials used to extract the data, we can see that 
the convergence is approximately fourth-order. 



tive impact on the numerical error. This, too, allows us 
to recover the expected (second-order) convergence. Fig- 
ures [S] and [HI illustrate the effects of the measures taken 
to control the convergence behavior. 
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FIG. 8: Behavior of convergence tests for a particle in circular 
orbit at r — 6 M. We show differences between simulations of 
the I — 2, m = 2 multipole moment using different step sizes 
of 2 (Via), 4 (</>4), 8 (ips), 16 (</>i 6 ), 32 (V> 33 ) and 64 (V> 64 ) cells 
per M. Displayed are the real part of the rescaled differences 
S4-2 = (ip4 — 1P2) etc. of the field values at the position of the 
particle, defined as in Fig. [5] The values have been rescaled 
so that — for fourth order convergence — the curves should all 
coincide. The upper panel corresponds to a set of simulations 
where the particle traverses the cells away from their tips. 
The curves do not coincide perfectly with each other, seem- 
ingly indicating a failure of the convergence. The lower panel 
was obtained in a simulation where the particle was carefully 
positioned so as to pass through the tips of each cell it tra- 
verses. This set of simulations passes the convergence test 
more convincingly. 



C. High-^ behavior of the multipole coefficients 



Inspection of Eq. (|4. 1 1|) reveals that a plot of as 
a function of £ (for a selected value of t) should display 
a linear growth in £ for large £. Removing the term 
should produce a constant curve, removing the term 
(given that Cr^ — 0) should produce a curve that decays 
as £~ 2 , and finally, removing the -D(^) term should pro- 
duce a curve that decays as £~ 4 . It is a powerful test of 
the numerical methods to check whether these expecta- 
tions are borne out by the numerical data. Fig. [TO] plots 
the remainders as obtained from our numerical simula- 
tion, demonstrating the expected behavior. It displays, 
on a logarithmic scale, the absolute value of Re ^(^w, the 
real part of the (+) component of the self-force. The orbit 
is eccentric (p — 7.2, e = 0.5), and all components of the 
self-force require regularization. The first curve (in trian- 
gles) shows the unregularized multipole coefficients that 
increase linearly in £, as confirmed by fitting a straight 
line to the data. The second curve (in squares) shows par- 
tially regularized coefficients, obtained after the removal 
of (£ + l/2)A( M }j this clearly approaches a constant for 
large values of £. The curve made up of diamonds shows 
the behavior after removal of B^y because Ci^ = 0, it 
decays behavior that is confirmed by a fit to 
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FIG. 9: Behavior of convergence tests for a particle in circular 
orbit at r = 6 M. We show differences between simulations 
of the I = 2, m = 2 multipole moment using different step 
sizes of 8 (tps), 16 (i/'ie), 32 (1/132), and 64 (064) cells per M. 
Displayed are the real part of the rescaled differences 5i6-8 = 
tpie ~ ^8 etc. of the field values at the position of the particle, 
defined as in Fig. [5] The values have been rescaled so that — 
for second order convergence — the curves should all coincide. 
The upper two panels correspond to simulations where the 
second order algorithm was used throughout. For the topmost 
one, care was taken to ensure that the particle passes through 
the tip of each cell it traverses, while in the middle one no 
such precaution was taken. Clearly the curves in the middle 
panel do not coincide with each other, indicating a failure 
of the second-order convergence of the code. The lower panel 
was obtained in a simulation using the mixed-order algorithm 
described in the text. While the curves still do not coincide 
precisely, the observed behavior is much closer to the expected 
one than for the purely second order algorithm. 



the £ > 5 part of the curve. Finally, after removal of 
D M /[(£ -\){£+ §)] the terms of the sum decrease in 
magnitude as £~ 4 for large values of £, as derived in [l5[. 
Each one of the last two curves would result in a con- 
verging sum, but the convergence is much faster after 
subtracting the Dr^ terms. We thereby gain more than 
2 orders of magnitude in the accuracy of the estimated 
sum. 

Figure [TO] provides a sensitive test of the implemen- 
tation of both the numerical and analytical parts of the 
calculation. Small mistakes in either one will cause the 
difference in Eq. (|4.11[) to have a vastly different behav- 
ior. 
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TABLE I: Results for the self-force on a scalar particle with 
scalar charge q on a circular orbit at ro = QM. The 
first column lists the results as calculated in this work us- 
ing time-domain numerical methods, while the second and 
third columns list the results as calculated in 0, HU using 
frequency-domain methods. For the t and 4> components the 
number of digits is limited by numerical roundoff error. For 
the r component the number of digits is limited by the trun- 
cation error of the sum of multipole coefficients. 





This work: 


Previous work: 


Diaz-Rivera 




time-domain 


frequency-domain [4] 


et. al. [14] 


All<j> R 

q t 


3.60339 x 10" 4 


3.60907254 x 10~ 4 




q r 


1.6767 x 10~ 4 


1.67730 x 10~ 4 


1.6772834 x 10~ 4 


q 


-5.30424 x 10~ 3 


-5.30423170 x 10~ 3 





FIG. 10: Multipole coefficients of the dimensionless self-force 



Re for a particle on an eccentric orbit (p — 7.2, e — 
0.5). The coefficients are extracted at t = 500 M along the 
trajectory shown in Fig. [12] The plots show several stages of 
the regularization procedure, with a closer description of the 
curves to be found in the text. 




is more efficient, and we expect the results of [1, [l4| to 
be much more accurate than our own results. This fact 
is reflected in the number of regularization coefficients 
we can reliably extract from the numerical data, before 
being limited by the accuracy of the numerical method: 
the frequency-domain calculation found usable multipole 
coefficients up to £ = 20, whereas our data for is 
dominated by noise by the time £ reaches 16. Figure [Til 
shows this behavior. 



E. Accuracy of the numerical method 

Several figures of merit can be used to estimate the 
accuracy of numerical values for the self-force. 

An estimate for the truncation error arising from cut- 
ting short the summation in Eq. 1|4.11[) at some £ max can 
be calculated by considering the behavior of the remain- 
ing terms for large £. Detweiler et. al. [HI showed that 
the remaining terms scale as £~ 4 for large £. They find 
the functional form of the terms to be 



FIG. 11: Multipole coefficients of $ (0) for a particle on a circu- 



lar orbit. Note that $ 



(0) 

is linked to $f via $f 



The multipole coefficients decay exponentially with 
I k, 16, at which point numerical errors start to dominate 



D. Self- force on a circular orbit 



I until 



EV- 



3/2 



(2£-3)(2^- l)(2£ + 3)(2^ + 5)' 



(5.3) 



where V3/2 = 36\/2. We fit a function of this form to the 
tail end of a plot of the multipole coefficients to find the 
coefficient E in Eq. (|5.3[) . Extrapolating to £ — > oo we 
find that the truncation error is 



For the case of a circular orbit, the regularization pa- 
rameters ^4(o) j B(o)j an( i -O(o) all vanish identically, so 
that the (0) (or alternatively the t) component of the 
self- force does not require regularization. Figure [TT| thus 
shows only one curve, with the magnitude of the multi- 
pole coefficients decaying exponentially with increasing 
I. 

As a final test, in Table UJ we compare our result for the 
self-force on a particle in a circular orbit at r = 6M to 
those obtained in 0, [T3 | using a frequency-domain code. 
For a circular orbit, a calculation in the frequency domain 



(5.4) 



12\/2£4 



(2£„ 



•3)(24 



l){2£n 



1)(24 



3) 

(5.5) 



where ^ m ax is the value at which we cut the summation 
short. For all but the special case of the (0) component 
for a circular orbit, for which all regularization parame- 
ters vanish identically, we use this approach to calculate 
an estimate for the truncation error. 
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A second source of error lies in the numerical calcula- 
tion of the retarded solution to the wave equation. This 
error depends on the step size h used to evolve the field 
forward in time. For a numerical scheme of a given con- 
vergence order, we can estimate this discretization error 
by extrapolating the differences of simulations using dif- 
ferent step sizes down to h = 0. This is what was done 
in the graphs shown in Sec. IV Bl 

We display results for mildly eccentric orbits. A high 
eccentricity causes (3 r <£> (displayed in Fig. [7} to be plagued 
by high frequency noise produced by effects similar to 
those described in Sec. IV Bl This makes it impossible to 
reliably estimate the discretization error for these orbits. 
We do not expect this to be very different from the errors 
for mildly eccentric orbits. 

Finally we compare our final results for the self-force 
F a to "reference values" . For circular orbits, frequency- 
domain calculations are much more accurate than our 
time-domain computations. We thus compare our results 
to the results obtained in Q. Tabic HT1 lists typical values 
for the various errors listed above. 



error estimation 


mildly eccentric orbit 


truncation error (-^-$(+)) 


« 2 x 1(T 3 % 


discretization error (^-d r &em) 


w 10" 5 % 


comparison with reference values 


circular orbit 


M' Z 77. 


0.2% 


q 2. 


0.04% 




2 x 10~ 4 % 



TABLE II: Estimated values for the various errors in the com- 
ponents of the self-force as described in the text. We show 
the truncation and discretization errors for a mildly eccentric 
orbit and the total error for a circular orbit. The truncation 
error is calculated using a plot similar to the one shown in 
Fig. 1161 The discretization error is estimated using a plot 
similar to that in Fig. [7] for the £ = 2, m = 2 mode, and the 
total error is estimated as the difference between our values 
and those of 0]. We use p = 7.2 , e = 0.5 for the mildly 
eccentric orbit. Note that we use the tetrad component $(+) 
for the truncation error and the vector component 9 r <E> for 
the discretization error. Both are related by the translation 
table Eqs. (|A6|) - (|A9|) . we expect corresponding errors to be 
comparable for $(+) and d r <&- 



VI. SAMPLE RESULTS 

In this section we describe some results of our numer- 
ical calculation. 

A. Mildly eccentric orbit 

We choose a particle on an eccentric orbit with p = 7.2, 
e = 0.5 which starts at r = pM/(l — e 2 ), halfway between 
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FIG. 12: Trajectory of a particle with p = 7.2, e = 0.5. The 
cross-hair indicates the point where the data for Fig. 1101 was 
extracted. 
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FIG. 13: Regularized dimensionless self-force ^-FV, -^4-F r 

and on a particle on an eccentric orbit with p = 7.2, 

e = q 5. 



periastron and apastron. The field is evolved for 1000 M 
with a resolution of 16 grid points per M, both in the t 
and r* directions, for £ — 0. Higher values of £ (and thus 
m) require a corresponding increase in the number of 
grid points used to achieve the same fractional accuracy. 
Multipole coefficients for < £ < 15 are calculated and 
used to reconstruct the regularized self-force F a along 
the geodesic. Figure H3] shows the result of the calcula- 
tion. For the choice of parameters used to calculate the 
force shown in Fig. 1131 the error bars corresponding to 
the truncation error (which are already much larger than 
than the discretization error) would be of the order of 
the line thickness and have not been drawn. 

Already for this small eccentricity, we see that the self- 
force is most important when the particle is closest to the 
black hole (ie. for 200 M <t< 400 M and 600 M <t< 
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FIG. 14: Trajectory of a particle on a zoom-whirl orbit with 
p = 7.8001, e = 0.9. The cross-hairs indicate the positions 
where the data shown in Fig. I16l and ll7l was extracted. 



800 M); the self-force acting on the particle is very small 
once the particle has moved away to r » 15 M. 



B. Zoom-whirl orbit 
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FIG. 15: Self-force acting on a particle. Shown is the dimen- 
sionless self- force ^4- Ft, ^-F r and HFa, on a zoom- whirl 
orbit with p — 7.8001, e = 0.9. The inset shows a magni- 
fied view of the self-force when the particle is about to enter 
the whirl phase. No error bars showing an estimate error are 
shown, as the errors shown eg. in Table [TT] are to small to 
show up on the graph. Notice that the self-force is essentially 
zero during the zoom phase 500 M < t < 2000 M and reaches 
a constant value very quickly after the particle enters into the 
whirl phase. 



Highly eccentric orbits are of most interest as sources 
of gravitational radiation. For nearly parabolic orbits 
with e < 1 and p > 6 + 2e, a particle revolves around the 
black hole a number of times, moving on a nearly circu- 
lar trajectory close to the event horizon ("whirl phase"), 
before moving away from the black hole ( "zoom phase" ) . 
During the whirl phase the particle is in the strong field 
region of the black hole, emitting copious amounts of 
radiation. Figures [14] and [15] show the trajectory of a 
particle and the force on such an orbit with p = 7.8001, 
e = 0.9. Even more so than for the mildly eccentric 
orbit discussed in Sec. IVI A[ the self-force (and thus the 
amount of radiation produced) is much larger while the 
particle is close to the black hole than when it zooms out. 

Defining energy E per unit mass and angular momen- 
tum L per unit mass in the usual way, 

and following eg. the treatment of Wald [16l |, Ap- 
pendix C, it is easy to see that the rates of change E 
and L (per unit proper time) are directly related to com- 
ponents of the acceleration a a (and therefore force) ex- 
perienced by the particle via 

E = —a t , L = a0. (6.2) 

The self-force shown in Fig. [15] therefore confirms our 
naive expectation that the self- force should decrease both 
the energy and angular momentum of the particle as ra- 
diation is emitted. 



It is instructive to have a closer look at the force acting 
on the particle when it is within the zoom phase, and also 
when it is moving around the black hole on the nearly cir- 
cular orbit of the whirl phase. In Fig. [TBI and Fig. [17] we 
show plots of $(o)£ vs. I after the removal of the Ar^, 
BffA, and D terms. While the particle is still zooming 
in toward the black hole, $>(o)e behaves exactly as for the 
mildly eccentric orbit described in Sec. IVI Al over the full 
range of £ plotted; ie. the magnitude of each term scales 
as £°, £~ 2 and £~ 4 , after removal of the At^, B^\, and 
Di^ terms respectively. Close to the black hole, on the 
other hand, the particle moves along a nearly circular tra- 
jectory. If the orbit were perfectly circular for all times, 
ie. r = 0, then the (0) component would not require reg- 
ularization at all, and the multipolc coefficients would 
decay exponentially, resulting in a straight line on the 
semi-logarithmic plot shown in Fig. [17] As the real orbit 
is not precisely circular, curves eventually deviate from a 
straight line. Removal of the A^ term is required almost 
immediately (beginning with £ ss 3) , while the term 
starts to become important only after I w 11. This shows 
that there is a smooth transition from the self-force on a 
circular orbit, which does not require regularization for 
the t and <fi components, to that of a generic orbit, for 
which all components of the self-force require regulariza- 
tion. 
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FIG. 16: Multipole coefficients of — Re &f s for a particle on 
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients are 
extracted at t = 2000 M as the particle is about to enter the 
whirl phase. As r is non-zero, all components of the self-force 
require regularization and we see that the dependence of the 
multipole coefficients on £ is as predicted by Eq. 11.91 After the 
removal of the regularization parameters v4( M ), -B( M ), and 
the remainder is proportional to £°, £~ 2 and £~ 4 respectively. 



APPENDIX A: TRANSLATION TABLES 

We quote the results of [4( for the translation table be- 
tween the modes and the tetrad components $( A1 )^ m 
with respect to the pseudo-Cartesian basis 
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, — cos 9 sin 
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r sin 
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r sin( 



and the complex combinations e?, •> := e?U ± ie 



'(2)' 



'(±) 



0, J] sinfle^, i cosfle* 4 * 
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With these, the spherical- harmonic modes $( A1 )£ m (t, r) 
are given in terms of <&^ m (i,r) by 
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FIG. 17: Multipole coefficients of Re for a particle on 
a zoom-whirl orbit (p = 7.8001, e = 0.9). The coefficients 
are extracted at t — 2150 M while the particle is in the whirl 
phase. The orbit is nearly circular at this time, causing the 
dependence on £ after removal of the regularization parame- 
ters to approximate that of a true circular orbit. 
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APPENDIX B: REGULARIZATION 
PARAMETERS 

For completeness we list the regularization parameters 
as calculated in Quantities bearing a subscript "0" 
are evaluated at the particle's position. 
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sign(A), 
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A, s — —p i 't l 0_ 
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where fo := 1 — 2M/tq and sign(A) is equal to +1 if 
A > and to -1 if A < 0. We have, in addition, A 
A {1) = Re[^ (+) ], and A {2) = lm[A {+) \. 
We also use 
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And finally, D ( _) = D {+) , D {1) = 
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£>(+) sin 0o - 
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APPENDIX C: PIECEWISE POLYNOMIALS 



(u + e, « ) 




(u , v + e) 



(to,r5) = (uo,i>o) 



(tt - e,w ) 



In two places in the numerical simulation we introduce 
piecewise polynomials to approximate the scalar field tpi m 
across the world line, where it is continuous but not dif- 
fer entiable . By a piecewise polynomial we mean a poly- 
nomial of the form 



( N 



p(t,r*) = { 



y Sip-u n v m iir*(u,v)>r* 



, (CI) 



y if r *(u,v) < r* 

^— ' n\nv 



where u — t—r*, v = t+r* are characteristic coordinates, 
Tq is the position of the particle at the time t(u, v), and 
N is the order of the polynomial, which for our purposes 
is N = 4 or less. The two sets of coefficients c nm and 
c' nm are not independent of each other, but are linked 
via jump conditions that can be derived from the wave 
equation [Eq. (|1.12[) ]. To do so, we rewrite the wave 
equation in the characteristic coordinates u and v and 
reintroduce the integral over the world line on the right- 
hand side, 

-4d u d v ip - Vip = [ S(T)8(u-u p )5(v-v p )dT, (C2) 



where S(t) — —8irq 



Y tm (ir/2,<l> (t)) 
Mr) 



is the source term and 



quantities bearing a subscript p are evaluated on the 
world line at proper time r. 

Here and in the following we use the notation 

[dZ8?1>] = Km [WV(*o, r* + e) - W^(*o, r* - e)] 

e— ►0+ 

(C3) 

to denote the jump in d™d™ip across the world line. First, 
we notice that the source term does not contain any 
derivatives of the Dirac J-function, causing the solution 
ip to be continuous. This means that the zeroth-order 
jump vanishes: [ip] = 0. Our task is then to find the re- 
maining jump conditions at a point (to, 7g) for n,m < 4. 
Alternatively, instead of crossing the world line along a 
line t = to = const we can also choose to cross along 
lines of u = uq = const or v = vq = const, noting that 
for a line of constant v the coordinate u runs from uq + e 
to uo — e to cross from the left to the right of the world 
line. Figure IT51 provides a clearer description of the paths 
taken. 



FIG. 18: Paths taken in the calculation of the jump condi- 
tions. (uo,vo) denotes an arbitrary but fixed point along the 
world line 7. The wave equation is integrated along the lines 
of constant u or v indicated in the sketch. Note that in order 
to move from the domain on the left to the domain on the 
right, u has to run from uo + e to uo — e. Where appropriate 
we label quantities connected to the domain on the left by a 
subscript "— " and quantities connected to the domain on the 
right by "+". 



In order to find the jump [d u ip] we integrate the wave 
equation along the line u = u from u - e to u + e 

-4 / d u d v ^dv - I Vipdv = 

J vq— e Jvq—c 

fVQ+€ 

S(t)S(uq — u p ) / 5(v — v p )dv dr, 

'7 J Do — e 



(C4) 



which, after involving J v °_ e S(v — v p )dv = 9{v p — vq + 
e)6(v -v p + e) and S(g(x)) = S(x - x )/ \g'(x )\, yields 

1 /o 



4E-f Q 



—S(t ), 



(C5) 



where the overdot denotes differentiation with respect to 
proper time r. 

Similarly, after first taking a derivative of the wave 
equation with respect to v and integrating from uq + e to 
Uq — e, we obtain 



110— e ru —e 

4 / d u dfydu - / 

luQ+t JUQ+t 

/•Uo— e 

S( T ) / S(u — u p )du5'(vo — v p )dr. 

'7 Juo+t 

(C6) 



We find 



rn2 ,1 1 fo d 
^^lEThd^ 



fp 



E + r 



-S(t) 



(C7) 



Systematically repeating this procedure we find expres- 
sions for the jumps in all the derivatives that are purely 
in the u or v direction. Table [TtTI lists these results. Jump 
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[fit] 



|t=t 



|t=t 



|t=t 



M = 



TABLE III: Jump conditions for the derivatives purely in the 
it or v directions, r and r are the particle's radial velocity and 
acceleration, respectively. They are obtained from the equa- 
tion of motion for the particle. £ := ^j 1 - and £ := ^j 1 - were 
introduced for notational convenience. Quantities bearing a 
subscript p are evaluated on the particle's world line, while 
quantities bearing a subscript are evaluated at the parti- 
cle's current position. Derivatives of V with respect to either 
u or v are evaluated as d u V = —\jd r V and d v V = \fd r V, 
respectively. 



conditions for derivatives involving both u and v are ob- 
tained directly from the wave equation [Eq. (|C2|) ]. We 
see that 

[8udvi>] = 0, (C8) 
and taking an additional derivative with respect to u on 



both sides reveals that 



(C9) 



Systematically repeating this procedure we can find jump 
conditions for each of the mixed derivatives by evaluating 

= -\ WWW)] , (cio) 

where n, m > and derivatives of V with respect to 
either u or v are evaluated as d u V = — \fd r V and d v V — 
hfd r V, respectively. 

The results of Table [m] and Eq. fCTO]) allow us to 
express the coefficients of the left-hand polynomial in 
Eq. (|C1|) in terms of the jump conditions and the co- 
efficients of the right-hand side: 



(Cll) 



For N = 4 this leaves us with 25 unknown coefficients 
c-nm which can be uniquely determined by demanding 
that the polynomial match the value of the field on the 
25 grid points surrounding the particle. When we are 
interested in integrating the polynomial, as in the case of 
the potential term in the fourth-order algorithm, we do 
not need all these terms. Instead, in order to calculate 
e.g. the integral ff y,Vtp dudv up to terms of order h 5 , 
as is needed to achieve overall 0(h 4 ) convergence, it is 
sufficient to include only terms such that n + m < 2, thus 
reducing the number of unknown coefficients to 6. In this 
case Eq. (|C1[) becomes 



Y ^HL u n v m iir*{u,v)>r* 
z — ' n'.m' 



E 

m+?i<2 



p(t,r*)=i m -^ 2 d . (C12) 



'm . 



mm 



The six coefficients can then be determined by matching 
the polynomial to the field values at the six grid points 
which lie within the past light cone of the grid point 
whose field value we want to calculate. 
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